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OPTIMIZATION  OF  CHEMICAL  VAPOR  INFILTRATION  WITH  SIMULTANEOUS 

POWDER  FORMATION* * * § 


A.  DITKOWSKI+,  D.  GOTTLIEB* ,  AND  B.W.  SHELDON^ 

Abstract.  A  key  difficulty  in  isothermal,  isobaric  chemical  vapor  infiltration  is  the  long  processing 
times  that  are  typically  required.  With  this  in  mind,  it  is  important  to  minimize  infiltration  times.  This 
optimization  problem  is  addressed  here,  using  a  relatively  simple  model  for  dilute  gases.  The  results  provide 
useful  asymptotic  expressions  for  the  minimum  time  and  corresponding  conditions.  These  approximations 
are  quantitatively  accurate  for  most  cases  of  interest,  where  relatively  uniform  infiltration  is  required.  They 
also  provide  useful  quantitative  insight  in  cases  where  less  uniformity  is  required.  The  effects  of  homogeneous 
nucleation  were  also  investigated.  This  does  not  affect  the  governing  equations  for  infiltration  of  a  porous 
body,  however,  powder  formation  can  restrict  the  range  of  permissible  infiltration  conditions.  This  was 
analyzed  for  the  case  of  carbon  infiltration  from  methane. 

Key  words,  composites,  chemical  vapor  deposition  (CVD),  optimization,  computer  simulation,  theory 

Subject  classification.  Physical  Sciences 

1.  Introduction.  A  variety  of  materials  are  produced  by  infiltration  processes.  In  these  techniques  a 
fluid  phase  (i.e.,  a  gas  or  a  liquid)  is  transported  into  a  porous  structure,  where  it  then  reacts  to  form  a  solid 
product.  These  methods  are  particularly  important  for  producing  composite  materials,  where  the  initial 
porous  perform  is  composed  of  the  reinforcement  phase  (i.e.,  fibers,  whiskers,  or  particles)  and  infiltration 
produces  the  matrix.  [1],  [2].  A  detailed  assessment  of  the  relevant  reaction  and  mass  transport  rates  during 
infiltration  requires  mathematical  modeling,  using  a  minimum  of  two  coupled  partial  differential  equations 
which  describe  changes  in  the  reactant  concentration  and  the  solid  structure  as  a  function  of  both  position 
and  time.  This  type  of  modeling  can  also  be  extended  to  analyze  the  optimization  and  control  of  infiltration 
processes. 

The  research  presented  here  specifically  considers  optimization  for  a  set  of  two  equations  which  describe 
isothermal,  isobaric  chemical  vapor  infiltration  (CVI).  In  this  process  a  vapor-phase  precursor  is  transported 
into  the  porous  preform,  and  a  combination  of  gas  and  surface  reactions  leads  to  the  deposition  of  the  solid 
matrix  phase.  During  infiltration  the  formation  of  the  solid  product  phase  eventually  closes  off  porosity  at 
the  external  surface  of  the  body,  blocking  the  flow  of  reactants  and  effectively  ending  the  process.  This  is  a 
key  feature  of  most  infiltration  processes.  Isothermal,  isobaric  CVI  often  requires  extremely  long  times,  so 
it  is  generally  important  to  minimize  the  total  processing  times. 

This  paper  considers  the  problem  of  determining  the  optimal  pressure  and  temperature  which  correspond 
to  the  minimum  infiltration  time.  From  a  practical  perspective,  the  nature  of  the  porous  preform  is  often 
predetermined  by  the  intended  application  (e.g.,  the  physical  dimensions  and  the  fiber  size  are  invariants). 
Thus,  the  process  can  only  be  controlled  with  process  variables:  temperature,  pressure,  and  gas  composition. 
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Note  that  the  pressure  and  temperature  determine  several  different  physical  quantities  in  the  model,  such  that 
the  general  understanding  of  the  optimum  conditions  is  not  immediately  obvious  from  the  basic  formulation. 

Several  previous  efforts  have  addressed  minimal  infiltration  times  [7],  [8],  [9],  [10],  [11].  The  work 
described  here  differs  from  these  previous  efforts  in  several  ways.  First,  the  asymptotic  solutions  obtained 
here  make  it  possible  to  estimate  the  optimal  conditions  without  conducting  detailed  numerical  calculations. 
This  paper  uses  a  relatively  simple  mathematical  model  which  is  based  on  the  diffusion  and  reaction  of 
a  single,  dilute  precursor.  However,  the  numerical  results  which  are  used  to  verify  and  understand  the 
limitations  of  the  asymptotic  results  are  generally  similar  to  previously  reported  numerical  results.  Also,  the 
results  presented  here  are  based  on  carbon  CVI  from  methane,  in  comparison  with  previous  efforts  on  CVI 
optimization  which  have  emphasized  SiC  CVI  from  methyltrichlorosilane.  The  use  of  different  chemistry  does 
not  change  the  general  conclusions  obtained  here,  however,  the  actual  values  do  correspond  to  a  somewhat 
different  process.  This  paper  also  considers  an  additional  constraint  due  to  homogeneous  nucleation  (powder 
formation),  which  has  not  been  treated  in  previous  efforts.  Powder  formation  can  be  significant  in  many 
CVI  processes  (including  carbon  formation),  and  this  effect  can  alter  the  optimal  temperature  and  pressure 
under  certain  conditions. 

This  paper  is  organized  as  follows:  Section  2  presents  the  basic  set  of  two  partial  differential  equa¬ 
tions  used  to  model  isothermal,  isobaric  CVI  (including  initial  and  boundary  conditions).  A  definition  for 
a  successful  process  and  a  discussion  on  the  optimization  problem  is  given.  In  Section  3  an  analysis  of  the 
optimization  problem  is  given.  The  analysis  performed  is  based  on  asymptotic  expansions  as  well  as  compu¬ 
tations.  The  results  of  the  analysis  are  optimal  working  pressure  and  temperature.  In  Section  4  the  effects 
of  powder  formation  in  the  analysis  are  included  in  the  analysis;  here  too  the  pressure  -and  temperature  to 
minimize  the  final  time  are  provided.  In  Section  5  a  discussion  of  the  significance  of  these  results  is  presented. 

2.  Formulation.  A  mathematical  description  of  infiltration  requires  one  or  more'  partial  differential 
equations  which  describe  the  evolution  of  the  matrix  (i.e.,  the  solid  phase),  and  one  additional  partial 
differential  equation  for  each  chemical  species  in  the  fluid  phase.  For  a  simple  pore  structure,  the  continuity 
equation  for  species  i  is 


(2.1)  ^  +  V.Ni  =  jTvirRr 

r 

where  t  is  time,  s  is  the  void  fraction  of  the  media,  C{  nd  N\  are  the  concentration  and  the  flux  of  species  i, 
nr  is  the  number  of  the  gaseous  species,  is  the  stoichiometric  coefficients  for  the  ith  gaseous  species  in 
the  rth  reaction,  and  Rr  represents  the  volumetric  reaction  rate  of  reaction  r. 

The  basic  partial  differential  equation  (s)  which  describe  reaction  and  mass-transport  in  porous  media 
(i.e.,  the  fluid  phase)  are  well-established  [3],  [5].  For  example,  the  Dusty-Gas  model  [6]  is  typically  used  to 
describe  multicomponent  diffusion  and  convection  in  a  porous  body. 

^  N<  .  RT  v  GjNi  -  CiNj  _  C{Be 

The  simplest  formulation  for  the  fluid  phase  is  obtained  by  considering  one  reacting  species.  For  highly 
diluted  reactant  systems,  the  Dusty-Gas  model  can  be  simplified  to  give  the  following  approximate  expression 
for  the  flux  in  one  spatial  dimension  ( Z ): 


(2.3) 


N  =  —D 


ac 

dZ 
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where  C  is  the  concentration  of  the  diluted  species  and  Z  is  the  distance  into  the  preform.  It  is  convenient  to 
write  C  in  terms  of  the  temperature  T,  the  total  pressure,  P,  and  the  mole  fraction  of  the  reacting  species, 
X: 


(2.4) 


C 


XP 
RT * 


For  a  diffusion-limited  process,  with  one  dilute  reactant  species  in  one  spatial  dimension,  Eq.  (2.1) 
becomes  (using  Eq.  (2.3)  and  Eq.  (2.4)  ): 


(2.5) 


d(sXP) 

dt 


d 

dZ 
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where  VM  is  the  molar  volume  of  the  solid  product  ,  u  is  the  rate  at  which  the  solid  product  grows  (vol¬ 
ume/area/time)  and  Sv(e)  is  the  gas/solid  surface  area  per  unit  volume  of  the  porous  solid.  The  last  term  in 
Eq.  (2.5)  describes  the  rate  at  which  the  gas-phase  precursor  is  consumed  (or  created)  by  chemical  reactions 
inside  of  the  pores,  with  the  assumption  that  there  are  no  homogeneous  gas-phase  reactions. 

Describing  the  evolution  of  the  matrix  phase  is  equivalent  to  considering  the  change  in  the  void  fraction, 
e  (i.e.,  the  volume  fraction  of  gas  inside  of  the  porous  solid).  The  evolution  of  e  is  given  by: 


(2.6)  g  =  -uSv(e). 

The  boundary  conditions  most  often  used  for  CVI  models  are  to  fix  the  concentration  at  the  outer  surface 
of  the  preform,  and  to  assume  that  diffusion  occurs  in  from  two  opposite  sides,  such  that  there  is  no  net  flux 
in  the  middle  of  the  preform  (i.e.,  at  Z  =  L,  where  L  is  the  half-thickness  of  the  preform): 


(2.7)  X(0  ,t)=X„ 

(2.8)  Xz(l,  t)  =  0 
where  z  =  Z/L . 

The  initial  condition  is  given  by: 


(2.9)  e(z,  0)  -e0{z). 

A  specific  CVI  model  requires  expressions  for  u ,  SVi  and  D.  Our  objective  in  this  work  is  to  use  relatively 
simple  formulations  for  each,  as  a  basis  for  assessing  optimization  during  isothermal,  isobaric  CVI.  As  an 
example,  consider  the  formation  of  carbon  matrix  composites  from  a  mixture  of  CH4  in  an  H2  carrier  gas, 
where  the  following  net  reaction  occurs: 


(2.10)  CH4  (g)  ^  C(s)  +  2  H2(g). 

The  form  of  Eq.  (2.5)  is  based  on  the  assumption  that  the  CH4  concentration,  C,  is  dilute  (i.e.,  the  reactant 
concentration  is  much  smaller  than  the  carrier  gas  concentration).  If  the  carbon  growth  rate  is  proportional 
to  the  precursor  concentration,  then: 


3 


(2.11) 


u 


where  k  is  the  reaction  rate  constant  (with  units  ~^)  and  R  is  the  gas  constant.  The  standard  Arrhenius 
expression  for  k  is: 


(2.12)  k  =  Akexp(-Q/T) 

where  Q  is  the  activation  energy  divided  by  the  R  and  Ak  is  a  pre-exponential  factor. 

The  preforms  used  for  CVI  typically  have  a  complex  porous  structure.  However,  a  cylindrical  pore  is 
often  used  to  formulate  simple  models.  This  leads  to  the  following  expression  for  Sv: 

(2.13)  Sv(e)  =  hEhU 

ro 

where  r0  is  the  initial  pore  radius  and  e0  is  the  initial  concentration  of  e. 

The  effective  diffusivity  of  the  dilute  species,  D,  can  be  expressed  as 

(2.14)  D=^DM[\  +  Nk(e)}-1 


Dm  — 


AmT 3/2 
P 


where  Dm  is  the  binary  diffusion  coefficient  for  the  reactant  species  in  the  carrier  gas  (e.g.,  CH 4  in  H2),  Am 
is  a  species  dependent  constant,  Nk  is  the  ratio  of  Dm  and  the  Knudsen  diffusion  coefficient,  and  6  is  the 
tortuosity  factor.  For  randomly  distributed  cylindrical  pores,  0  is  estimated  to  be  3,  and  Nk  is  given  by: 


where 


Nk 


AdT 

rjlf2P 


(2.15) 


_  AmtjqT^ 

D%P 


D°k  is  the  Knudsen  diffusion  coefficient  for  the  initial  pore  size,  at  some  reference  temperature  T*. 
Substituting  Eqs.  (2.11)  and  (2.13)  into  Eqs.  (2.6)  and  (2.5)  gives  the  following  forms: 


drjjzji) 


dt 


d_ 

dz 


dc(z,t)\ 
dz  ) 


-fc(z,t) 
a2  7j(z,t)c(z,t) 


(2.16) 

(2.17) 
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Table  2.1 

Values  of  the  constants. 


■Aq, 

6.35  (10)21  K3/2  atm"1 

A0 

3.64  (10)18  K3/2  atm"1  sec"1 

Vo  = 

o 

II 

•+o 

0.85 

Q 

55000  K  [17] 

Ad' 

10/xm  diam  fibers 

1.54  (10)-5  atm/K 

80 nm  diam  fibers 

1.93(10)-®  atm/K 

200Mm  diam  fibers 

7.7  (10) -7  atm/K 

where 

(2.18) 

t!  \  P?73 

f  M  ~  3(Pt7  +  AdT) 

(2.19) 

77  =  v/e 

(2.20) 

,  pe(~Q/T) 

a  ~ Aa  7,3/2 

(2.21) 

X0Pe(-°/T) 

P  =  - 

(2.22) 

_  2 AkJT0L2 
a  Vm^oAmR 

(2.23) 

A  2Aky/eo 

0  ~  r0R 

Where  7*0  is  the  initial  pore  radius  and  Xo  denotes  the  fraction  of  the  active  gas  in  the  inlet  (at  z  —  0).  Note 
that  a 2  is  dimensionless  and  that  /?  has  units  of  inverse  time.  The  time  derivative  in  Eq.  (2.5)  has  been 
removed  in  Eq.  (2.17),  which  is  permissible  because  a  pseudo-steady  state  C  profile  is  achieved  in  a  short 
amount  of  time  (i.e.,  compared  to  the  time  scale  over  which  e  changes)  [4],  Transforming  e  to  77  simplifies 
equation  (2.16).  Note  that  77  is  proportional  to  Sv,  so  rj  can  be  viewed  as  a  dimensionless  surface  area  per 
volume.  Values  of  the  constants  in  Eqs.(2.18)-(2.23)  are  given  in  Table  2.1,  for  the  case  of  carbon  CVI  from 
methane,  (see  Eq.(2.10)),  with  a  preform  thickness  of  2L  =  3mm. 

The  system  (2.16),  (2.17)  is  subject  to  the  initial  condition 


(2.24) 


*7(2,0) =770(2). 


This  paper  treats  only  the  case  of  a  uniform  initial  condition,  770(2)  =  770.  The  boundary  conditions  are  (see 
(2.7),  (2.8)): 


(2.25) 

(2.26) 


c(0,  t)  =  1 
dc 


dz 


=  0. 


Z=1 


With  this  model  it  can  be  shown  that: 

•  The  value  of  the  void  function  in  the  inlet  z  —  0,  is 


(2.27) 


t?(0,i)  -  ~t 
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•  There  exists  a  critical  time  tc  =  ^p-.  At  this  time,  the  void  function  vanishes  at  z  =  0,  the  inlet 
closes  completely,  and  the  process  ends. 

•  The  void  function  rj  and  the  concentration  function  c  are  bounded  from  above  and  below.  0  < 
c(z,t)  <  1,  r}(Q,t)  <  r)(z,t)  <  rjo  for  t  <tc  —  ^p. 

•  For  any  time  t  <  tCi  the  void  function  r}(z,t)  is  monotonically  increasing  function  of  the  spatial 
variable  2. 

•  The  concentration  function  c(z,t)  is  monotonically  decreasing  function  of  the  spatial  variable  2. 

The  literature  on  diffusion  and  reaction  in  porous  media  typically  uses  a  dimensionless  ratio  of  the  diffu¬ 
sion  and  reaction  rates,  sometimes  known  as  the  Thiele  modulus.  This  parameter  varies  during  infiltration, 
because  of  changes  in  the  microstructure  with  time.  Thus,  previous  CVI  modeling  has  used  an  initial  Thiele 
modulus  as  an  approximate  assessment  of  the  relative  infiltration  kinetics  [4],  [7],  [8].  In  terms  of  the  for¬ 
mulation  specified  here,  the  initial  Thiele  modulus  is  equal  to  a2??(2,0)/ 7(77(2, 0)).  In  general,  when  a2  is 
small,  diffusion  is  fast  and  infiltration  is  relatively  uniform.  When  a2  is  large,  the  deposition  reaction  is  fast 
and  infiltration  is  highly  non-uniform. 

The  parameters  a2  and  /?  depend  on  the  three  key  process  variables:  T,  P,  and  Xq.  Process  optimiza¬ 
tion  during  CVI  is  achieved  by  setting  these  variables  to  optimal  values.  In  isothermal,  isobaric  CVI  the 
infiltration  kinetics  are  controlled  by  diffusion  and  the  deposition  reaction.  To  achieve  relatively  uniform 
infiltration,  diffusion  must  be  fast  relative  to  the  deposition  rate.  This  is  typically  accomplished  by  choos¬ 
ing  processing  conditions  that  result  in  a  slow  deposition  rate,  which  usually  leads  to  long  infiltration  times. 
Thus,  the  primary  basis  for  process  optimization  is  to  obtain  the  desired  amount  of  infiltration  in  the  shortest 
possible  time. 

A  general  definition  of  a  successful  process  includes  two  considerations: 

1.  At  the  end  of  the  process  (  at  time  t  —  tf)  the  void  function  77(2,  tf)  should  be  a  small  fraction  of 
its  initial  value,  either  in  the  whole  interval  or  in  a  certain  portion  of  the  interval  0  <  2  <  z\. 

2.  For  the  process  to  yield  good  results  it  is  important  that  the  void  function  is  uniformly  small  along 
the  2  axis. 

Mathematically,  we  express  these  considerations  by  stating  that  a  process  is  successful  if  for  some  time  tf 

(2.28)  7]{zx ,tf)  <  hr] 0 

(2.29)  r)(Q,tf)  =  k0r)o 

ko  «  1,  ko  <  k\.  Equation  (2.28)  states  that  the  final  values  of  the  void  function  77  should  be  small  in 
the  interval  between  the  inlet  and  the  point  z\  (note  that  77(2,  t)  is  monotonically  increasing  function  of  the 
spatial  coordinate  2).  In  most  problems  of  interest  21  —  1.  Conditions  (2.28)  and  (2.29)  state  that  the  void 
function  should  be  uniformly  small.  Also,  from  (2.27),  (2.29),  the  final  time,  £/,  is  given  by: 

(2.30)  tf  =  (l~ko)m^ 

Note  that  the  time  for  the  process  to  end  decreases  as  a  function  of  f3  (itself  a  function  of  the  temperature 
and  pressure,  given  in  (2.21)). 

The  goal  of  the  analysis  in  the  following  sections  is  to  find  the  temperature  and  pressure  that  minimize 
the  final  time  tf  for  achieving  a  successful  process. 
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3.  Analysis  of  the  Optimization  Problem.  Evaluation  of  Eq.  (2.27)  requires  an  approximation  for 
r)(z\,tf)  in  terms  of  the  pressure  and  temperature.  There  are  two  ways  to  obtain  it:  either  by  asymptotic 
expansions  or  by  direct  numerical  solution  of  the  system  (2. 16), (2. 17). 


3.1.  Asymptotic  Expansions.  The  set  of  equations  (2.16),  (2.17)  is  nonlinear  and  explicit  solutions 
can  not  be  obtained.  However  some  analytical  approximation  may  be  derived  by  noting  that  in  most  of  the 
problems  of  interest  a2  is  small,  and  therefore  it  makes  sense  to  expand  the  solutions  in  power  of  a2.  The 
details  of  obtaining  this  expansion  to  order  a2  will  be  given  elsewhere  [16].  For  small  a2  the  final  result  is: 


(3.1) 

(3.2) 
Where 

(3.3) 


c(zyt)  ~  1  -  a2(z  -  z2/ 2) 


E 


r}(Zyt)  ~  E(t)  4-  3 a2(z  -  z2/ 2) 


m = 


It  can  be  shown  that  the  error  in  the  expansion  is  proportional  to  j  ,  where  fa  is  the  ratio  between  the 
final  void  function  at  the  inlet  and  its  initial  value,  as  defined  in  (2.29).  In  the  next  subsection  the  numerical 
results  and  the  asymptotic  expansions  are  compared,  to  demonstrate  the  validity  of  the  expansions  in  the 
range  of  relevant  a2. 

The  asymptotic  expansion  (3.2)  is  used  to  get  an  explicit  form  for  the  uniformity  constraint  (2.28)  in 
terms  of  the  temperature  T  and  the  pressure  P.  Substituting  (3.2)  into  (2.28)  gives: 


(3.4) 


Mo  +  3a2  (zi  -  zl/2) 


<  Mo 


Substituting  a2  from  (2.20)  and  rearranging  gives: 


(3.5) 


J(P,T)  <  (fa  -  fa)r)0. 


where 

106  + 

Rearranging  terms  shows  that  uniformity  is  assured  if 


(3.6)  J(P,T)  =  3Aa^e-^T(Zl  -  z2J 2) 


(3.7) 
where 

(3.8) 

(3.9) 

(3.10) 


p  <  ^0T3/2 eQ/T  _ 
“  #2 


Bo 

b2 


T, 


B0  = 


Bi 


rjojh  -  kp) 
AaZ(zi  -  z\l 2) 

1) 


Ap .  1 
Vo 


These  results  can  now  be  used  to  approximate  the  temperature  and  pressure  that  minimize  the  infiltration 
time.  Recall  (see  (2.30))  that  the  final  time  tf  is  inversely  proportional  to  /3  given  in  (2.21).  The  final  time 
tf  is  therefore  minimized  if  the  function 


7 


(3.11) 


F(P,T)  =  -e~Q/T 

is  maximized  under  the  uniformity  constraint  (2.29).  Inspection  of  (3.11)  shows  that  in  order  to  maximize 
we  have  to  take  the  equality  sign  in  (3.7).  Substituting  this  into  (3.11)  ,  it  is  easily  verified  that  the  following 
function  must  be  maximized: 


(3.12) 


G(T)  =  ^T1/2  -  ^ e~Q/T 
B2  -D2 


This  indicates  that  the  final  time  tf  to  achieve  a  successful  process  is  minimized  by  choosing  temperature 
and  pressure  satisfying 


(3.13) 

(3.14) 


T3/2eQ/T  _ 

Bo 


P=|i(Q-T). 


(3.15) 


Moreover  the  minimal  final  time  t™171  is  given  by 

2(1  -  k0)  \og(l/k0)  AaQ 


±mm  _ 

tf  ~ 


hi  —  ko  ApX  o 

Where  the  temperature  T  is  given  by  (3.13). 


3(*i-*?/2) 


1 


TV2(Q-T) 


The  explicit  formulas  (3.13)-(3.15)  lead  to  the  following  observations: 

1.  The  minimum  final  time,  t7Jlzn)  decreases  as  Ap  decreases.  (For  example,  as  the  molecular  diffusion 
becomes  dominant.) 

2.  t™in  decreases  as  k\  increases.  This  reflects  the  fact  that  increasing  k\  relaxes  the  uniformity 
condition. 

3.  As  z\  increases  toward  z\  —  1,  the  minimum  final  time  t™m  increases. 

3.2.  Computational  Results.  In  the  previous  subsection  the  asymptotic  expansion  of  rj  was  used  to 
define  a  functional  J(P,T)  such  that  each  pair  P,  T  that  satisfies  J(P,  T)  <  ( ki  —  ko)  rjo  leads  to  a  solution 
that  satisfies  the  conditions  for  a  successful  process,  (2.28),  (2.29).  The  optimal  P  and  T  was  then  obtained 
based  on  the  final  time.  This  result  is  approximately  correct  since  the  asymptotic  expansion  was  used  to 
approximate  condition  (2.28).  This  section  uses  numerical  solutions  of  (2.16),  (2.17)  to  create  a  ‘numerical 
J  functional5,  (i.e.,  a  functional  relation  between  P  and  T  that  ensures  a  successful  process). 

Two  algorithms  were  used  to  solve  this  problem:  one  based  on  a  finite  difference  approximation  and  one 
on  spectral  methods.  These  are  described  in  the  Appendix.  The  schemes  were  run  with  fco  =  0.1,  z±  =  1  and 
ki  =  0.15  or  0.7.  Note  that  ki  =  .15  corresponds  to  relatively  uniform  infiltration  while  k\  =  .7  is  relatively 
non-uniform.  Although  most  applications  require  relatively  uniform  infiltration  (i.e.,  lower  fci),  there  are 
some  cases  where  a  non-uniform  profile  may  be  desirable.  Two  reasons  for  a  higher  k\  are  that  it  enables 
faster  infiltration  times,  and  it  produces  materials  with  lower  density.  For  example,  both  of  these  attributes 
are  desirable  during  the  formation  of  thin  carbon-carbon  composites  for  bipolar  plates  in  proton  exchange 
membrane  (PEM)  fuel  cells  [18]. 

Three  values  were  taken  for  A/),  1.54  (10)-5,  1.93  (10)“6  and  7.7  (10)“7  (see  Table  2.1).  Plots  of  the 
numerical  and  the  asymptotic  J  curves  are  presented  in  Figure  3.1,  plots  of  tf  vs.  P  are  presented  in 
Figure  3.2. 
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Fig.  3.1.  J -Curves  obtained  both  numerically  (solid  lines)  and  by  asymptotic  analysis  (dotted  lines).  Conditions  on  or 
below  these  lines  will  satisfy  the  uniformity  condition  (Eq.  (3.35)).  Conditions  above  this  line  do  not  satisfy  the  uniformity 
condition.  All  cases  here  correspond  to  Xq  —  0.1,  with:  (a)  k\  =  0.15,  10/im  diameter  fibers;  (b)  k\  =  0.15,  200/um  diameter 
fibers;  (c)  ki  =  0.7,  lOjim  diameter  fibers;  (d)  k\  =  0.7,  200/im  diameter  fibers. 


Pressure  (atm)  Pressure  (atm) 


Fig.  3.2.  Times  which  correspond  to  conditions  on  the  J -curves  in  Figure  3.1,  with  values  obtained  both  numerically 
(solid  lines)  and  by  asymptotic  analysis  (dotted  lines).  The  filled  circles  show  the  minimal  time,  (a)  k\  —  0.15,  10^m  diameter 
fibers;  (b)  ki  —  0.15,  200/im  diameter  fibers;  (c)  k\  =  0.7,  10/xm  diameter  fibers;  (d)  fci  =  0.7,  200^m  diameter  fibers. 

The  process  is  very  sensitive  to  changes  in  the  temperature  T,  as  can  be  seen,  for  example,  by  comparing 
the  results  for  P  =  .latm,T  =  1387 K,  vs.  P  =  .latm,  T  =  1448 K  (for  k\  =  0.7  and  Ad  —  1.54  (10)-5). 
Increasing  temperature  by  61  degrees,  (~  4%),  decreases  the  final  time  from  160.3min  to  31.4min,  i.e.,  by  a 
factor  of  5  and  produces  an  infiltration  profile  which  is  much  less  uniform.  This  occurs  because  of  the  strong 
temperature  dependence  of  the  deposition  reaction.  As  the  fibers  size  increases  tf  decreases  slightly  and  the 
minimal  time  occurs  at  lower  pressures. 

When  the  uniformity  requirement  dictates  that  k\  —  .15  in  (2.28),  the  condition  J(P,T)  —  (&i  -  ko)rjo 
yields  a  P  and  T  pair  such  that  a2  ~  .01.  In  this  case,  the  asymptotic  expansions  agree  well  with  the 
numerical  results.  The  predicted  temperature  that  assures  uniformity  differs  by  only  few  degrees  from  the 
one  obtained  numerically,  and  tf  differs  by  less  than  10%.  When  ki  was  increased  to  0.7,  The  large  value  of 
a2 ,  (^  0.1)  leads  to  more  inaccuracy  in  tf,  (See  Table  3.1).  However,  the  optimal  temperatures  predicted 
by  the  asymptotic  result  is  still  accurate  within  10%. 
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Table  3.1 

t f  for  the  optimal  P  and  T. 


M D 

h 

P  [atm] 

Numerics 

T[K] 

tf  [H]  | 

Asymptotics 
P  [atm]  T[K] 

tf  [H] 

1.54  (lO)"5 

0.15 

6.605 

1201.65 

16.03 

7.694 

1195.06 

17.62 

1.93  (10) -6 

0.15 

0.830 

1258.48 

15.27 

0.962 

1251.37 

16.80 

7.7  (lO)-7 

0.15 

0.334 

1389.75 

14.93 

0.373 

1282.43 

16.42 

1.54  (lO)"5 

0.7 

3.185 

1338.99 

0.3934 

7.694 

1265.75 

1.428 

1.93  (lO)"6 

0.7 

0.375 

1492.76 

0.3727 

0.962 

1405.28 

1.357 

7.7  (10)  "7 

0.7 

0.155 

1573.00 

0.3632 

0.373 

1478.51 

1.324 

In  all  cases,  the  asymptotic  results  agree  qualitatively  with  the  numerical  results.  The  curves  obtained 
numerically  were  almost  parallel  to  the  asymptotic,  and  the  points  of  minima  are  almost  in  the  same  place. 
The  asymptotic  results  are  conservative,  they  always  overestimated  the  final  time,  and  gave  more  restrictive 
conditions  on  P  and  T  for  uniformity,  i.e.,  P  and  T  obtained  by  the  asymptotic  analysis  never  predict  a 
successful  process  if  it  does  not  exist.  However,  since  P  and  T  obtained  by  the  asymptotic  analysis  may  be 
much  different  from  the  optimal  ones  (obtained  by  the  numerical  analysis),  tf  may  be  much  larger  then  the 
minimal  value  (See  Table  3.1). 

4.  Homogeneous  Nucleation.  CVI  processes  can  be  limited  by  homogeneous  nucleation  (i.e.,  powder 
formation)  in  the  gas  phase.  This  effect  has  not  been  treated  in  previous  CVI  models  because  it  generally 
occurs  outside  of  the  solid  preform.  However,  powder  formation  can  impose  serious  limitations  on  CVI 
operating  conditions  during  the  formation  of  carbon  and  oxide  matrices.  Thus,  this  phenomena  imposes  a 
constraint  on  the  allowable  CVI  operating  conditions.  Nucleation  kinetics  are  typically  described  with: 


(4.1)  I  =  ZrrA'iXoP/RT)  exp(-A g*/kBT) 

where  I  is  the  steady- state  nucleation  rate,  A*  is  the  area  of  a  critical  cluster,  J*  is  the  flux  at  which 
atoms  are  added  to  a  critical  cluster,  A g*  is  the  free  energy  barrier  to  forming  a  critical  cluster,  ks  is 
Boltzmann’s  constant,  and  Zj  is  the  so-called  Zeldovich  factor,  see  [12].  A  rigorous  model  requires  that  Zj 
be  evaluated  numerically,  however,  using  standard  approximations  for  Z/,  the  pre-exponential  terms  in  (4.1) 
can  be  combined  to  give: 


(4.2)  I  Sf  (X0P/RT)22(y/kBT)^2Ak  exp (-Q/RT)  exp(-Ag*/kT) 

where  7  is  the  surface  free  energy  of  the  cluster,  and  the  constants  Ak  and  Q  describe  the  reaction  rate 
constant  (see  (2.12)).  In  practice,  the  permissible  value  of  I  depends  on  the  reactor  configuration,  as  well  as 
its  actual  value.  For  the  current  analysis,  we  assume  that  powder  formation  limits  CVI  when  the  nucleation 
rate  exceeds  some  allowable  level,  Ium.  With  this  in  mind,  Eq.  (4.1)  can  be  revised  to  yield: 
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Fig.  4.1.  Effect  of  homogeneous  nucleation  for  k\  =  0.7,  Xq  =0.01,  10/xm  diameter  fibers,  (a)  The  numerically  obtained 
J -Curve  (uniformity  constraint,  Eq.  (3.5))  and  the  I-curve  (nucleation  limit,  Eq.  (4-3)  with  A]\m  =  2.0  1021/  (b)  Limiting 
time  as  a  function  of  pressure.  The  left  part  of  the  curve  is  determined  by  the  numerically  obtained  J -curve  and  the  right  side 
is  determined  by  the  I-curve,  with  the  minimal  time  shown  by  the  filled  circle.  The  dotted  line  corresponds  to  the  approximate 
J -curve  which  was  determined  with  asymptotics. 


(4.3)  Aum  =  himk'ff/Ak  7  =  (X0P)2/r2'5  exp  (-Aj/T™) 

where  the  two  exponential  terms  in  Eq.  (4.2)  have  been  combined  to  give  one  term  with  two  empirical 
constants,  Aj  and  m .  Since  A g*  can  have  a  relatively  complex  T  dependence,  this  empirical  approach 
was  adopted  to  provide  a  relatively  simple  expression.  This  form  was  applied  to  the  results  of  Loll  et 
al.,  who  report  threshold  conditions  for  the  onset  of  significant  nucleation  (i.e.,  Xo  vs.  T  at  P  —  latm) 
[13],  [14].  A  good  fit  to  their  experimental  data  was  obtained  with  m  —  1.5,  Aj  =  750,000  Km,  and 
Aum  —  3.3  (10) ”17  atm2/K2-5.  In  general,  the  value  of  Aum  is  somewhat  arbitrary,  since  it  reflects  a 
threshold  for  a  given  reactor.  By  varying  Aum ,  it  is  possible  to  assess  different  tolerance  levels  for  powder 
formation.  For  example,  recent  carbon  CVI  experiments  at  Oak  Ridge  National  Laboratory  tolerate  higher 
powder  formation  levels  than  those  described  by  Loll  et  al.  with  a  threshold  value  that  corresponds  to 
Aum  =  2.0  (10)-21  atm2/K2,5  [15]. 

The  effect  of  adding  the  powder  formation  constraint  can  be  seen  from  Figure  4.1,  where  the  /-curves 
are  defined  by  Eq.  (4.2).  As  seen  in  Figure  4.1a,  the  new  constraint  limits  the  pressures  and  temperatures 
to  values  which  are  below  both  the  I  and  the  J-curves.  For  a  given  pressure,  the  minimal  tf  corresponds 
to  a  temperature  on  either  the  I  or  J- curve  (whichever  is  lower) .  Thus  if  the  minimal  tf  found  in  Section  3 
(i.e.,  when  P  and  T  are  on  the  J-curves)  to  the  left  of  the  J- curve,  then  the  additional  constraint  does  not 
change  the  previous  results.  If,  on  the  other  hand,  this  point  is  on  the  right  of  the  I- curve  then,  the  minimal 
tf  occurs  at  the  intersection  between  the  I  and  J-curves.  This  point  can  be  clearly  seen  as  a  cusp  in  the  tf 
vs.  P  in  the  right  plot. 

A  complete  assessment  of  X0  effects  requires  solutions  with  the  full  Dusty-Gas  model,  because  large 
values  of  Xo  violate  the  assumption  of  a  dilute  reactant  gas.  However,  considering  only  values  up  to  Xo  =  0.1 
provides  useful  insight  into  optimizing  dilute  systems.  Without  the  homogeneous  nucleation  constraint  (i.e., 
as  Aum  —►  oo),  the  minimal  time  is  inversely  proportional  to  Xo,  and  the  optimal  pressure  and  temperature 
do  not  vary  with  Xq  (see  Section  3).  However,  homogeneous  nucleation  limits  the  operating  conditions  when 
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Fig.  4.2.  Effect  ofX o  with  10pm  diameter  fibers:  (a)  optimal  pressures;  (b)  optimal  temperatures ;  (c)  minimum  infiltration 
times  and  with  200pm  diameter  fibers:  (d)  optimal  pressures ;  (e)  optimal  temperatures ;  (f)  minimum  infiltration  times. 


Aum  is  low  enough,  as  illustrated  in  Fig  3.  The  effect  of  this  limitation  on  the  optimal  conditions  and  on 
the  minimal  time  are  shown  in  Figs.  4  and  5.  These  results  lead  to  the  following  conclusions: 

1.  As  in  Section  3,  the  asymptotic  results  are  in  good  agreement  with  the  numerical  results  for  ki  =  .15, 
where  a2  ~  .01,  and  much  less  accurate  in  the  case  k\  —  .7  where  a2  ~  .1.  In  both  cases,  however, 
the  asymptotic  results  agree  qualitatively  with  the  numerical  results. 

2.  Notwithstanding  the  dilute  reactant  gas  restriction,  tf  are  monotonically  decreasing  functions  of 
Xo,  thus  it  is  advisable  to  work  in  the  ‘highest’  Xo  possible.  However  since  the  optimal  P  is  also 
a  monotonically  decreasing  functions  of  Xo,  this  value  of  X0  is  limited  by  the  lowest  operational 


12 


200  A 


Fig.  4.3.  Effect  of  Alim  on  tf,  with  k\  =  0.15;  (a)  10pm  diameter  fibers  and  Xo  —  0.1;  (b)  10/um  diameter  fibers  with 
Xq  —  0.02;  (c)  10j/m  diameter  fibers  with  Xo  ~  0.1.  All  values  are  based  on  numerical  results. 

a)  b) 


Fig.  4.4.  Effect  of  Ayim  on  the  optimal  pressure  (a)  and  temperature  (b)}  for  the  same  cases  plotted  in  Fig.  5. 


pressure.  For  example  for  working  pressure  of  about  0.01atm,  the  maximum  allowable  Xo  is  only 
0.05,  for  200Mm  diameter  fibers  (Ad  ~  7.7  (10)”7  atmK-1 ). 

3.  For  a  given  Xq,  as  the  fiber  diameter  increases,  P  and  tf  decrease  and  T  increases.  But  unlike 
Section  3,  the  differences  here  are  significant.  This  occurs  because  the  homogeneous  nucleation 
condition  forces  us  to  work  in  a  region  where  the  dependence  on  Ad  is  much  stronger. 

4.  The  homogeneous  nucleation  constraint  causes  the  optimal  temperature  and  pressure  to  vary  with 

X0. 

As  the  value  of  Aum  increases,  the  effects  of  homogeneous  nucleation  are  less  severe.  This  can  be  seen 
in  Fig.  6,  which  shows  the  effects  of  varying  Ium.  Note  that  the  minimum  infiltration  time  is  dramatically 
increased  when  there  is  a  significant  limitation  imposed  by  homogeneous  nucleation.  In  general,  the  process 
must  be  operated  at  lower  pressures  to  avoid  powder  formation.  Some  increase  in  the  corresponding  optimal 
temperature  accompanies  this  decrease  in  pressure.  The  slope  discontinuities  in  Fig.  6  correspond  to  the 
conditions  where  the  homogeneous  nucleation  constraint  no  longer  has  an  effect,  (i.e.,  the  optimal  conditions 
are  determined  solely  by  the  J  curve). 
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5.  Conclusions.  Minimizing  infiltration  times  for  isothermal,  isobaric  CVI  is  important  because  pro¬ 
cessing  times  are  typically  very  long.  The  work  presented  here  provides  a  detailed  assessment  of  the  pressure 
and  temperature  which  will  minimize  the  total  required  time,  based  on  a  simplified  model  for  a  single,  dilute 
reactant  species.  This  formulation  makes  it  possible  to  understand  the  basic  physics  of  the  problem  in  terms 
of  a  relatively  small  number  of  lumped  parameters.  The  basic  objective  of  this  optimization  problem  is  to 
obtain  a  density  profile  with  a  prescribed  amount  of  uniformity,  in  the  shortest  possible  time  (Section  3). 
The  asymptotic  results  are  particularly  useful,  because  they  make  it  possible  to  determine  optimal  conditions 
without  doing  numerical  calculations  (under  conditions  where  a2  is  small  enough).  Based  on  comparisons 
with  the  numerical  results,  the  asymptotic  forms  are  also  qualitatively  accurate  when  a 2  is  larger.  Thus, 
the  asymptotic  results  provide  a  clear  understanding  of  how  the  optimal  conditions  are  related  to  the  key 
parameters  for  the  problem. 

The  effects  of  homogeneous  nucleation  were  also  analyzed,  as  an  additional  constraint  on  the  basic 
optimization  problem.  This  issue  has  not  been  considered  in  previous  work  on  CVI  modeling,  however, 
it  can  limit  operating  conditions  in  systems  were  powder  formation  is  significant  (e.g.,  the  formation  of 
carbon  matrix  composites).  The  results  obtained  here  provide  a  quantitative  assessment  of  the  conditions 
where  homogeneous  nucleation  imposes  limitations  on  infiltration  conditions.  When  these  limitations  occur, 
powder  formation  also  increases  the  minimum  infiltration  time. 


6.  Appendix.  To  gain  confidence  in  our  computations,  two  completely  different  numerical  methods 
were  used:  the  pseudospectral  method  and  a  finite  difference  method. 

•  In  the  pseudospectral  Chebyshev  method  the  grid  points  is  chosen  to  be 


(6.1) 


0  <  j  <  N 


where  N  is  the  total  number  of  grid  points  and  was  10  for  most  of  the  runnings. 

The  spectral  differentiation  matrix  takes  the  value  of  a  given  function  at  the  g  rid  points  Zj  and 
yields  the  values  of  the  derivative  of  the  interpolation  polynomial  at  these  points.  The  points  Zj 
are  the  nodes  of  the  Gauss  Lobatto  Chebyshev  quadrature  formula.  The  matrix  can  be  written 
explicitly: 


Djk  - 

1  C,  (-l)i+l 

j¥=k, 

2  c*  sin  27 jr(i+*)  sin 

(6.2) 

Djj  — 

1  Zj 

2  sin *($i) 

j  /  0,  AT 

Doo  = 

-dnn  =  2JV«+1 

We  apply  the  matrix  Djk  twice,  once  for  the  vector  c  taking  into  account  the  boundary  condition 
c(0,t)  =  1  and  then  to  f(rj)  cz  and  taking  into  account  that  cz(l,t)  —  0.  This  yields  a  linear  system 
for  the  values  of  c(zj,t). 

In  the  next  stage  we  update  rj  by  the  standard  third  order  Runge-Kutta  scheme. 

•  A  second-order  finite-difference  scheme  using  the  equidistance  grid 

(6.3)  Zj  =  L=jh  0  <  j  <  N  . 
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The  differentiation  matrix  can  be  written  explicitly: 


DX1 

—  —■fiz(f{zi—h/2)  +  f(zi  +  h/2)) 

D12 

=  T?f(zi+h/2) 

Djj-i 

II 

sh 

dT 

! 

>• 

Djj 

Djj+i 

=  +  V2) 

Dn  n~i 

cT 

1 

r-H 

II 

Dn  n 

=  -p-  (/(l  ~ /i/2)  +  »?) 

j  =  1 


j  =  2...N  —  1 


f  in  the  mid-points  was  interpolated.  In  each  step  we  solve  the  system 


f  2)  \ 


Dc  = 


In  the  next  stage  we  use  c  to  update  r)  by  the  standard  fourth  order  Runge-Kutta  scheme.  Since 
this  scheme  is  less  accurate  we  used  80  grid  points. 

The  results  of  both  schemes  were  compared  and  the  differences  in  r)(z,tf)  were  less  then  10-6. 
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